shockname = {'Housing supply','Housing','Mortgage','Expectation';
            'shock','consumption shock','rate shock','shock'};
shockname2 = {'Housing supply shock','Housing consumption shock', 'Mortgage rate shock','Expectation shock'};


             
colorcode=[0 0.5 0; 0 0 1; 0 1 1; 1 0 0;0.4 0.4 0.4;];

meany= repmat(mean(y,1),size(y,1),1);


 
[T,n_var] = size(y);

n_rent=n_var+1;
n_pr_cum=n_var+2;
n_gdp_cum=n_var+3;
%n_inv_r=n_var+4;
n_hh_form_cum=n_hh_form+4;
series{n_pr_cum}='Price Level';    
%series{n_inv_r}='Residential Investment Level';    
series{n_gdp_cum}='GDP Level';    
series{n_rent}='Rent';    
series{n_hh_form_cum}='Cum. net HH form';    

 %%%%%%%Set size of first Block (zero for no zero restrictions)
%n_chol=0;



%% Estitmation 
 [gamma,Su,resid,Par] = BVARgibbs_parameters(y,p,tau); %OLS Values;
 var_Ferr=diag(Su); % Variance of forecast errors

%% Identification
    tel=0; %% Counts number of accepted draws
    tel2=0;
    rot=0; % Number of draws
    nonstab_tot=0; % Number of unstable draws (at least one Eigenvalue >1)
        while tel<n_rot_acc   % Fix number of acceptance            
%for jj = 1:n_draws
jj=1;
    if paramuc==1;  % with parameter uncertainty (for Figure 2)   
    nonstab=1;
    while nonstab==1; % only stable draws
    [gamma_g(:,:,jj),Su_g(:,:,jj)]=BVARgibbs_draws(Par,1);
    [choleskyimp(:,:,:,jj),A(:,:,jj), nonstab] = IRF_BVARGibbs(gamma_g,Su_g,p,n_var,1,n_step); %This is sampling whith parameter uncertainty
    nonstab_tot=nonstab_tot+nonstab;
    end
    else
    [choleskyimp,A ] = IRF_BVARGibbs_2(gamma,Su,p,n_var,1,n_step); %% This is sampling whithout parameter uncertainty
    gamma_g=gamma;
    end
           
          %Draws Rotation matri x
            W = normrnd(0,1,n_var,n_var);
                                    Q= getqr(W);
      for tt=1:n_step
            irfcand(:,:,tt)=choleskyimp(:,:,tt,jj)*Q;  %jj refers to parameterdraw jj, tt is the step in the IRF
      end
      

 
irfcand(n_rent,:,:)=irfcand(n_rent_price,:,:) +cumsum(irfcand(n_price,:,:),3);
irfcand(n_pr_cum,:,:)=cumsum(irfcand(n_price,:,:),3);
irfcand(n_gdp_cum,:,:)=cumsum(irfcand(n_gdp,:,:),3);
irfcand(n_hh_form_cum,:,:)=cumsum(irfcand(n_hh_form,:,:),3);
J = zeros(n_var,n_var*p); J(:,1:n_var) = eye(n_var);


        for shock=1:n_var; 
            % (1) Supply shock
         
               test_su= [checksign3(irfcand(n_pr_cum ,shock,:),LL,KK,'+')  
                checksign3(irfcand(n_total_vac, shock,:),LL,KK,'-')   
                checksign3(irfcand(n_perm ,shock,:),LL,KK,'-')      
                ];
          
             if min(test_su)==1
                 shockcheck(shock,1)=1;
             elseif max(test_su)==-1;
             shockcheck(shock,1)=1;
             irfcand(:,shock,:)=-irfcand(:,shock,:);
             else
                shockcheck(shock,1)=0;
             end
            % (2) Demand Shock
         
              test_de=[  checksign3(irfcand(n_pr_cum, shock,:),LL,KK,'+')
                checksign3(irfcand(n_perm,shock,:),LL,KK,'+')
                checksign3(irfcand(n_total_vac, shock,:),LL,KK,'-')  
                checksign3(irfcand(n_mrate, shock,:),LL,KK,'+') 
                ];
            if min(test_de)==1
                 shockcheck(shock,2)=1;
             elseif max(test_de)==-1;
             shockcheck(shock,2)=1;
             irfcand(:,shock,:)=-irfcand(:,shock,:);
             else
                shockcheck(shock,2)=0;
            end
             
            % (3) Mortgage rate shock
             test_int=[
                checksign3(irfcand(n_pr_cum,shock,:),LL,KK,'+')
                checksign3(irfcand(n_perm ,shock,:),LL,KK,'+')     
                 checksign3(irfcand(n_mrate ,shock,:),LL,KK,'-')   
               ];
           
          if min(test_int)==1
                 shockcheck(shock,3)=1;
             elseif max(test_int)==-1;
             shockcheck(shock,3)=1;
             irfcand(:,shock,:)=-irfcand(:,shock,:);
             else
                shockcheck(shock,3)=0;
          end
            
          
            % (4) Expectation shock
         test_sp=[
                checksign3(irfcand(n_pr_cum,shock,:),LL,KK,'+')          
                checksign3(irfcand(n_perm ,shock,:),LL,KK,'+')     
                checksign3(irfcand(n_total_vac, shock,:),LL,KK,'+')  
                checksign3(irfcand(n_mrate ,shock,:),LL,KK,'+')      
                ];
            
            if min(test_sp)==1
                 shockcheck(shock,4)=1;
             elseif max(test_sp)==-1;
             shockcheck(shock,4)=1;
             irfcand(:,shock,:)=-irfcand(:,shock,:);
             else
                shockcheck(shock,4)=0;
            end
            
             
        end

        controlvec= sum(shockcheck,1);
    
            
        if controlvec==ones(1,size(controlvec,2));
             tel=tel+1;
             tel2=tel2+1;
             if tel/10==fix(tel/10) % Show number of accepted draws every 10 accepted draws
             disp([num2str(tel) ' accepted draws out of ' num2str(rot) ' draws' '('  num2str(nonstab_tot) 'unstable)']);
             end
            [irfstr(:,:,:,tel2), vardecomp(:,:,:,tel2)] = irfselec(shockcheck,irfcand);
            BB(:,:,tel2)=irfcand(1:n_var,1:n_var,1);
           errstr=BB(:,:,tel2)\resid';
            strerr(:,:,tel2) = ord_strerr(shockcheck, errstr);
            
        end

    
   clear irfcand
   rot=rot+1;
        end

rejectrate=100*(rot-tel)/rot; 
disp(['TOTAL DRAWS: ' num2str(rot) ' , TOTAL ACCEPTED: ' num2str(tel) ', REJECTION RATE: ' num2str(rejectrate)]);

no_str=size(shockcheck,2);                                       %Number of Structural Shocks 
    
%%Transformations
%Which Variables should be cumulated to convert from quarter to quarter growth to Level in IRF 
cum2=[n_price n_gdp];




%Transform into real investment level
irfstr_cum=irfstr;
% Transform into level to later convert to annualized growth
irfstr_cum(cum2,:,:,:)=cumsum(irfstr(cum2,:,:,:),3);
irfstr_cum(n_hh_form,:,:,:)=irfstr(n_hh_form_cum,:,:,:);

%
irfstr_mean=mean(irfstr_cum(:,1:no_str,:,:),4);
irfstr_cum_mean=mean(irfstr_cum(:,1:no_str,:,:),4);


irfstr_median=median(irfstr_cum(:,1:no_str,:,:),4);
irfstr_cum_median=median(irfstr_cum(:,1:no_str,:,:),4);

strerr_save=strerr;

